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Abstract. 

This  paper  presents  a  class  of  new  explicit  second  order  accurate  finite 
difference  schemes  for  the  ccnputaiton  of  weak  solutions  of  hyperbolic  conservation 
laws.  These  highly  nonlinear  schemes  are  obtained  by  applying  a  nonoscillatory 
first  order  accurate  schene  to  an  appropriately  modified  flux  function.  The  so 
derived  second  order  accurate  5?oh??Tes  achieve  high  resolution  while  preserving 
the  robus-tness  of  the  original  nonoscillatory  first  order  accurate  scheme. 

K\imerical  experiments  are  presented  to  demonstrate  the  performance  of  these 
new  schemes. 
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1.   Introduction. 

In  this  paper  we  consider  numerical  approximations  to  weak  solutions  of  the 
initial  value  problem  (IVP)  for  hyperbolic  systems  of  conservation  laws. 


u^  +  f(Li)j^  =  0,   uCx,0)  =  (|>(x),   -»  <  X  <  ".  (1.1) 


Here  u(x,t)  is  a  clourm  vector  of  m  unknowns,  and  f  (u) ,  the  flux,  is  a 
vector  valued  function  of  m  components.  Cl.l)  is  called  hyperbolic  if  all 
eigen'^talues   a-,(u),...,  a^Cu)  of  the  Jacobian  matrix  ACu) 


ACu)  =  f^  a.  2a) 


are  real  and  the  set  of  right  eignevectors  R^(u),...,  I^(u)  is  complete. 

Ke  assume  that  the  eigenvalues  {a  (u)}  are  arranged  in  a  nondecreacing  order 


a^Cu)  s:  a^(u)  <  ...  ^  a'^Cu) . 


(1.2b) 


We  consider  systems  of  conservation  laws  (1.1)  that  possess  an  entropy 
function  U(u),  defined  as  follows: 


(i)  U  is  a  convex  function  of  u,  i.e.,  U   >  0. 

'      '   uu    ' 


(ii)  U  satisfies 


vhere  F  is  sai^e  other  function  called  entropy  flux.  Admissible  weak  solutions 
of  (1.1)  satisfy,  in  the  weak  sense,  the  following  inequality: 
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U(u)^  +  FCu)^  ^0  a.3b) 


(see  [11]).  The  inequality  (1.3b)  is  called  an  entropy  condition. 

In  the  following  we  shall  discuss  numerical  approximations  to  weak  solutions 
of  Q.l)  vW.ch  are  obtained  by  C2k+  l)-point  explicit  schemes  in  conservation 
fomi 

3      D     [3+7    D-yJ 
where 

Vfere  V.  =  vCjAx,  nAt),  an  f  is  a  numerical  flux  function.  We  require  the 
numerical  flux  function  to  be  consistent  with  the  flux  f(u)  in  the  following 
sense: 

f(u,...,u)  =  f(u).  (l.Uc) 

Vfe  say  that  the  difference  scheme   (l.U)  is  consistent  with  the  entropy  condition 
(1.3b)   if  an.  inequality  of  the  following  kind  is  satisfied: 

lf}^h  if}  -  X    li^i     -  F?_l|  a.5a) 

viiere  U^  =  U(Vj),  7.^1  =  F(v"  ^^, , . .  ,v^_^^)  ;  here  ?  is  a  numerical  entropy 
flux,  consistent  with  the  entropy  flux  F(u),  i.e. 
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F(u,...,u)  =  FCu).  Cl.5b) 

Vfe  turn  now  to  discuss  the  question  of  convergence  of  the  finite  difference 
solution  of  Q.4)  to  weak  solutions  of  the  conservation  laws  (1.1).  Since  the 
finite-difference  scheme  is  nonlinear  and  the  conputed  solutions  are  certainly 
not  smooth,  therefore  L„-stability  of  a  consistent  finite-diference  scheme 
does  not  irriDly  convergence.  One  can  establish  convergence  of  finite  -  difference 
solutions  of  Cl.4)  to  weak  solutions  of  Q.l)  when  the  following  conditions 
are  satisfied: 

(i)  The  total  variation  with  respect  to  x  of  the  finite-difference  solutions 
is  uniformly  bounded  with  respect  to  t,  At  and  Ax. 

Cii)  The  finite-difference  scheme  (1.4)  is  consistent  with  the  entropy  condition 
Q.3b)  for  all  entropy  functions  of  Q.l). 

-Ciii)  The  entropy  condition  (1.3b)  implies  uniqueness  of  the  solution  to  the  IVP 
(1.1). 

Using  compactness  arguments  one  can  deduce  from  condition  (i)  the  existence 
of  convergent  s'ibsequences .  The  conservation  form  (1.4)  and  condition  (ii) 
imply  that  each  limit  solution  is  a  weak  solution  which  satisfies  the  entropy 
condition  Q.3b).  When  the  entropy  condition  inplies  uniqueness  of  the  IVP 
(conditon  (iii))  then  all  subsequences  have  the  same  limit  solution,  and  con- 
sequently the  finite-difference  scheme  is  convergent,  (see  [2]  ,  [9]  ,  [10]  )  . 

It  seems  possible  to  satisfy  conditions  (i)  and  (ii)  by  adding  a  hefty 
araount  of  ai"tificial  viscosity  to  the  finite  difference  scheme  (1.4).  The  addi- 
tional viscosity  terrr^s  darrp  possible  oscillations  in  the  computed  solution, and 


-  u  - 


HBke  the  ccfnvergence  process  simulate  the  zero-dissipation-limit  which  is  used 
to  select  the  unique  physically  relevant  weak  solution.  Unfortunately,  viscosity 
represents  an  irretrievahle  loss  of  inforriHtion  and  therefore  the  addition  of 
artificial  viscosity  brings  about  some  deterioration  in  resolution. 

In  this  paper  we  describe  a  new  nethod  to  design  finite- difference  schemes 
that  satisfy  conditions  Q)  and  Cii),  but  are  second-order  accurate  and  have 
high  resolution. 
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2.  V£ynotordcit--/   in  the  scalar  case. 

In  this  section  v;e  consider  the  IVP  for  a  scalar  conservation  law. 

"t  ■"  ^^""^x  -  ""t  ■"  ^^''^''x  "  °'   ^^""^  ^  aU^  ^^'^^^ 

u(x,0}   =  <.(x),  -CD  <  X  <  oo,  C2.1b) 

where    0 (x)    is   assumed   to   be   of   bounded   total   variation.       Every 
weak    solution   of    the    scalar    IVP    (2.1)    which   satisfies    the   entropy 
condition   has    the    following   mono  tonicity   property   as   a    function   of   t: 

i)       No  new  local  extrena  in     x     may  be  created. 

ii)     The  value  of  a  local  minimum  is  nondecreasing,  the  value  of  a  local  majdjnum 
is  nonincreasing . 

It  folla^7s  from  this  monotonicity  property  that  the  total  variation  in 
X,  T/(u(t)),    of     uCxjt),      is  nonincreasing  in  t,      i.e., 

1V(u(t2))    ^  TVCu(t^)),     for  all  t^  >.  t-^.  (2.2) 

Vfe  consider  now  explicit  (2k  +  1)  point  finite  difference  schemes  in  con- 
servation forrc  (1.4)  apprcxiiTHting  (2.1) 


n+1   TT,  n    n         n  ^    n  ,r  t^  n        ,n  . 

■j  --  ^^(Vj.k'  ^j-kn'---'^j+k)  =  ^3  -''^'^^^:-k.i'-'-'"j+k> 

(2.3a) 


^^3_5,-•-v.^^_^)] 


and  denote  (2.3a)  in  an  ooerator  form  as 
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v"^^  =  L-v""  C2.3b) 


We  say  that  the  finite  difference  scheme  (2.3)  is  total  variation  nonincreas- 
ing  (TV?JI)  if"  for  all  v  of  bounded  total  variation 

TV(L-v)  c   T/Cv)  C2.Ua) 

v.'here 

TVCu)  =  i_Uj+^/o^l      ;  (2.4b) 

here,  and  throughout  this  paper,  we  use  the  standard  notation 


^3+1/2^=  ^jn-^j'  ^2.5) 


Vfe  say  that  the  finite-difference  scheme  (2.3)  is  monotonicity  preserving 
^   thie  finite  difference  operator  L  is  monotonicity  preserving,  i.e.,  if  v 
is  a  monotone  miesh  function  so  is  L'V. 

We  say  that  the  finite  difference  scheme  (2.3)  is  a  monotone  scheme  if 
H  in  (2.3a)  is  a  monotone  nondecreasing  function  of  each  of  its  2k  +  1 
arguments. 

The  following  theorem  states  the  hierarchy  of  these  properties. 

Theoi^em.  2.1.   (i)  A  monotone  scheme  is  TVNI.   (ii)  A  IVNI  scheme  is  monotonicity 
pr>eser\'ing . 
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Proof:  (i)  It  was  proven  by  B.  Keyfitz  in  [8]  that  monotone  schemes  form  an 
t-, -contractive  semigroup,  i.e. 

IlL-v-  L-z|l   .<  I|v-  zll  (2.6a) 

CO 

for  all  i, -sunrable  v  and  z;  Here  I  Iu|  L  =  I   |u.j.   (2.4)  follows 

-L  ^1   j=-»  ^ 

imniediately  frcm  applying  (2.6a)  to  v  and  z  =  T-v,  (i.e.  z.  =  v.^^  for  all  j), 

(ii)  Let  (2.3)  be  a  TViJI  scheme  and  let  v  be  a  monotone  mesh  function  of 
bounded  total  variation,  and  denote  w  =  L*v.  Since  L  l-as  a  finite  support 
of  2k  +  1  points  it  is  sufficient  to  prove  that  w  is  nonotone  for  all  v  of 
the  form 

constant  =  v-,  Ii  <■  '^^ 

monotone  J_  ^  j  ^  J^     j  J+  ^  '^_ 

constant  =  v^  j  ^  J^. 

TV(v)  =  K  -  ^lI  ^2. 6b) 

V7e  prove  (ii)  by  negation.  Suppose  w  is  not  monotone,  then  it  has  at 

least  one  local  irixdinum  and  one  local,  maximum.  Denote  by  v  and  v^  the  values 

m      M 

of  the  first  two  successive  local  extrema,  then 

TV(w)  >.  |vj^  -  v^l  +  \v^  -v^\    >   TV(v), 

which  contn^dicts  the  assurrption  that  the  scheme  is  TVNI.     This  coirpletes  the 
proof  of  Tr.eorem     2.1 


hbnotone  scheiies  approximate  solutions  of  the  viscous  modified  equation 


u^  +  f(u)^  =  At[e(u,X)  u^]^   '  ^  =  ^  ^^-"^^^ 


e( 


(u,X)  =  ^-r  I     l\(iu,u,...,u)   -  X  a  Cu) 


2.^  ,    i. 


(2.7b) 


2A  f^^-k 
6(w,X)  >.   0,   5Cu,X)  I  0  C2.7c) 


to  second  order  accuracy;  Since  3(u  A)  |  0,  monotone  schemes  are  necessarily 
first  order  accurate;  H^   in  C2.7c)  denotes  |rr^w_,  ,  w  ,^  . . .  ,w,  ) .   (  see  [  8  ]) 


Since  monotone  schemes  are  T/>JT,  there  exist  convergent  subsequences  for  all 
initial  data  of  bounded  total  variation.  Each  lirrit  is  a  weak  solution  of  (2.1) 
that  satisfies  Oleinik's  entropy  condjition  (see  [  8  ]).  Since  OleiniJc's  entropy 
condition  inplies  uniqueness  of  the  IVP  (2.1),  we  conclude  that  all  subsequences 
converge  to  the  same  liTrit,  and  therefore  the  scheme  is  convergent  (see  [2  ]). 

Let  us  consider  new  the  scalar  constant  coefficient  case  a(u)  =  constant 
in  (2,1).  A  linear  finite-difference  approxiriHrion 

k 
v.   =  y  c.v.,„      ,  c.=  const.  (2.8a) 

is  monotonicity  preserving  if"  and  only  if 


c^  0  ,  -k  <:  £  <:  k. 
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(see  [4  ]).  Hence  any  linear  monotonicity  preserving  scheme,  and  therefore  any 
TT.n.     linear  schene,  is  a  monotone  scheme,  and  consequently  first  order  accurate. 

We  remark  that  the  previous  statement  does  not  exclude  the  possibility  of 
having  nonlinear  monotonicity  preserving  and  TVNl  schemes  that  are  second  order 
accurate  Cand  consequently  are  not  monotone  schemes).  In  fact  the  schemes  presented 
in  17  2  and  [6  J  are  monotonicity  preserving  (at  least  in  the  constant  coefficient 
case)  and  second  order  accurate. 

It  is  the  purpose  of  this  paper  to  present  new  high-resolution  second  order 
accurate  TVNI  scheires.  These  new  schemes  are  generated  by  converting  known 
3-point  first  order  accurate  TYfJI  schemes  into  new  5 -point  second  order  accurate 
TVNI  schemes.   Both  the  3-point  schemes  and  the  new  5-point  schemes  can  be  re- 
written ifi  the  form 


v"^l  =  L-v"  ^2. 9a) 


A.  ,  ,„v  (2.9b) 


(L-v).  =  Vj  .  C,^.,.^^,^^^/^v  -  C_^._^^2^._^^2 


where  A._^]^  is  defined  in  (2.5)  and 


^+,jn/2  =  C^(Vj_-Vj,v^+;^,v.^2^,   C_^^_^/2  =  C_(v._2,v._^,v.,v.^^)  (2.9c) 

Tne  follCTs-ing  Lemra  states  conditions  on  the  coefficients  (2.9c)  which  are 
sufficient  to  ensoure  that  the  schene  (2.9  )  is  TVNI. 
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Leama  2.2 

Let  the  coefficients  C      in   (2.9c)    satisfy  the  inequalities 

C_,j+i/2lO.        C^,j+i/2iO  (2.10a) 

^-.j+1/2  ■'■^+,3+1/2     ^^      '  (2.10b) 

then  the  scheme  (2.9)  is  TVNI. 

Proof:   Denote  w  =  L.v  and  subtract  (2.9b)  at  j  =  i  from  (2.9b)  at  j  =  i+1 
to  obtain 


^i+1/2"  =  C-.i-l/2'^i-l/2^  "■   ^^-^-,i+l/2-C+.i+l/2>^i+l/2^  +  ^+,1+3/2  ^i+3/2^- 


(2.11) 


By  (2.10)  all  the  coefficients  in  (j2.ll)  are  non  negative,  therefore 
l^i>l/2"'  ^'  ^1-C-, in/2  -^+,1+1/2^1^+1/2^'  "C-,i-l/2''i-l/2^1 


(2.12) 


*   C+,i+3/2l'i+3/2^' 


Sunning  C2.12)  for   -co  <  i  <  »  ^g  obtain 

*  i   =-  i-l/2l'i-l/2"l  *    J     <:.,i«/2l'i*3/2''l  =  J  J'in/2^1  '   ■"'^' 

vihlch   shows  that  (2.U)  is  satisfied.  The  equality  is  obtained  by  changing  the 
suiTYTHtion  index  in  the  last  t\-.'o  sums  in  the  RHS  of  tr;  inequality. 
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In  the  next  section  V7e  shall  use  heiTZTS.   2.2  to  design  second  order  accurate 
T/KL   sch^Ties.  Vfe  rer-ark  that  any  3  point  finite  difference  schaiie  in  conservation 
form  with  a  different iable  numerical  flux  can  be  rewritten  as  (2.9),  in  the  follow- 
ing way:  It  follows  from  the  ir^an  value  theorem  that  there  exist  C_^     and  C_ 
such  that 

Xlf(Vj,Vj,3_)  -  fCv.,v.)]  =  -C,(v.,v.,^)  A.,^/2^  (2.13a) 


X[f(v^_-L,Vj)  -  f(Vj,v^)]  =  -  C_(Vj_^,Vj)  Aj_^/2V   .  (2.13b) 


Expressing  the  numerical  flux  values  in  C2.3a)  vrLth  k  =  1  by  (2.13)  results 
in  the  form  C2.9). 


-  12  - 


3.       Second  order  accurate     TV!II     scheines. 


Let  us  consider  a  general  3-point     finite  difference  scheme  in  conservation 
form  (l.U)  with  a  numerical  flux     f    of  the  form 


fCv.,v.,^)  4lfCVj)  .  f(v.,^)  -  lQCAi.^|)   A.^1  V] 


v^ere 


^H 


[f(v.,,)-fCVj)]/A.,lv 


aCvj) 


when     A. ,1  v  4  0 
l'*'T 


when     A. ^1  V  =  0 
3+T 


CS.la) 


(3.1b) 


Here     Q(x)  is  some  function,  which  is  often  referred  to  as  the  coefficient  of 
numerical  viscosity. 

LenriH  3. 1     Let     Q(x)   in  (3.1a)   satisfy  the  fxiequalities 


|xl    $  Q(x)    <  1 


for  0   «:    |x|    ^  y   <:  1, 


(3.2) 


then  the  finite-difference  scheme  (1.4)  with  (3.1)  is  TVNI  under  the  CFL-like 
restriction. 


Xmax .  I  a'?  a  I  ^  P 
J  J  2 


(3.3) 


Proof:  Using  the  notation 


v.^1  -  .Xa.^l 
1  +  9     J+T 


(3.4a) 
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where     a.,]^     is     (3.1b),     v;e  rewrite   (3.1a)     as 
^2 


>-f.^l  =  U-(v.,v.^^)   =  Xf(v.)   -  7  -v.,1  .  Q(».,1V 


A.    1  V 


(3.Ub) 


and  similarly 


Xf.  1  =  >.f(v.  ,,v.)  =  Af(v.)  -  i 
:-2     D-1  :      3    ^ 


}  _1  +  0(v._l) 

J  9        J  9 


A.  1  V 


C3.Uc) 


Substituting  (3.4)  for  numerical  flux  values  in  (l.U)  we  get  the  form  (2.9) 


n+1    n 


n  ,  1| 


=  v"  -  X(f.  1  -  f .  1)  =  vj  +  yFcX"^.  ,1)  -  v.^ 


L   ^^2 


A._,l  V 
3+T 


(3.5a) 


WS-"-^  *  \] 


A.  1  v"  =  v"  +  C^,.,l  A.^1  v^  -  C  ,.  1  A.  1  v". 


2   -  2 


2   ^  2 


C    1  -  i 


Ip-Vi^*'i4. 


Since 


C.,._,l  +  C  ,.^1  =  Q(  V. 
+3+9    -  3+0      3+ 


1) 


(3.5b) 


(3.5c) 


it  follou's  from  C3.2)  and  (3.5)  that  corditions  (2.10)  of  Lemrra  2.2  are  satisfied 
under  the  CFL  restriction  (3.3)  and  therefore  the  finite-difference  schems  (3.1) 
is  TVNI. 

Tne  second  order  accjrate  Lax-Vfendpoff  scheme  has  the  num.erical  flux  (3.1) 


with  Qix)   =   X  ,   i.e. , 
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J.W 


Clearly  a  numerical  flux  of  a  second  order  accurate  scheme  f-;+^ /o  ^'^^^ 
to  satisfy 

fj+l/2  -  ^"1/2  =  «(''''  "-^^ 

for  all  siTooth  solutions  of  (1.1);  here  A  is  the  discretization  parameter. 
\Ihen  QCx)  is  constrained  by  C3.2),  then  the  3-point  scheme  (3.1)  is  only 
first  order  accurate,  for 

^3.1/2  =  ^"1/2  -  W  f-<^j.l/2>'  *   Q'^M/2"  'in/2  ""  ".7) 

and  therefore 

lfn/2  -  fj.l/2l    >■    h   '!^3*l/2l    -   '^jtl/2>'3-|'3+:/2  ^"l    =  °<"- 

Vte  describe  now  how  to  convert  a  S-point  first  crdsr  accurate  TVNI  scheme  to 


a  5-point  second  order  accurate  T'\/NT  schame.  Consider  the  application  of  a 
3-point  first  order  accurate  T\ 
of  the  original  flux  f Cu) :  Set 


^ 

3-point  first  order  accurate  TV!-!!  scheme  (3.1)  to  modiried  mesh  values  f^ 


^  =  f(v.)  +  ^  g,  ,  g,  =  g(Vi_3  ,v,.v,^,)  (3.8a) 

where  g  is  a  function  which  will  be  specified  below.  Tie  modified  numerical 
flux  f^^,/2  =  f^(v._^,v.,v._^j^,v._^2^  ^^  obtained  by  rcrlacing  f(v^),  ^(v^+j^) 
in  (3.1a),  (3.1b),  by  the  modified  values: 
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where 

^i+1/2   =   ^+1/2   -^  ^i+1/2    '  ^i+1/2   =    ^gl+l-8i^^'i+l/2  ^-  3-^"^ 

We  can   rewrite   (3.8b)    as 

^j+1/2   =i    [f(v^)+f(Vj^l)5    -^  YX    fSj+gj+i-Q(Vl/2^  Vl/2^    Vl/2  ^^-         ^'•''^^ 
Lpmma   3.2      Suppose  0(x)    is   Lipschitz   continuous  and   g        satisfies 


^^  (3.9b) 


T3n/2'^j+i/2  ^  ^gj.i  -  Sj  =  °^^  ^ 


then  the  nuirerical  flux  (3.8)   satisfies   C3.6) 


-M 

Proof:      The  rrodified  riLner-ical  fl-jx     f^+i/o      (-3.8d)     differs  from  the  original 
flux     f..      C3.1)     in  the  follot:ing  way  : 

(3.10a) 


Substit-oting  (3.7)  for  f.,-,/o  in  (3.10a)  we  get  that  (3.6)  holds  if  the 
^  3+1/2 

relation 

Cgj  ^  gj^3_)  -  [QCv^.i/o)  -  QCvj,i/2  ^-  -^3^/2^^    ^jn/2  ^    = 

(3.10b) 

f^Q^^3^1/2^   -    ^^3^/2^'^   ^3^/2  -  ^  °^^'^ 
is  satisfied. 
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Since     QCx)     is  Lipschitz  continuous 

"therefore  if  follows  from  (3.9b)  that  the  second  term  on  the  IJiS  of  (3.10b) 

2 

IS  itself  0(A  );  consequently  (3.9a)  implies  (3.10b).  This  conpletes  the  proof 

of  Leimra  3.2. 

We  construct  g.  =  g(v.  ,,  v., v.,,)  that  satisfies  (3.9)  in  the  following 
t.-j_    fe  x-1    i'  1+1  ^ 

way: 

g.  =  s.^^/2"^[0,mir.(|^.,^/2l^  gi-1/2 -^1+1/2)  ]  ^^"^^^^ 


'^i+l/2'"^^|gi+l/2l'    ISi-i/2i)  when     gi+i/2-gl-i/2   >'  ^ 

'V. 

i+1/2'%-1/2 


0  when     g,^,.„ •!.,,_   i  0 


where 


gi+1/2   =    7   tQ^^i+1/2^   -   ^Vl/2>'^  ^+1/2  -  ^3-1^^ 


^i+1/2   =   ^^^   *^i+l/2^-  ^^-^^^^ 

Lenma  3.  3:     Let     g.     be  defined  by  (3.11) , then  relations  (3.9a)  and  (3.9b)  are 
satisfied,  and 


l^j-M/21    =    1%-H  -  Sjl/l^jn/2-l    ^^    7I     ^5^^3+1/2^   -   ^^jn/2^'l-  ^3.12) 
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Proof:    First   let  us  assume  thia"^    '^^+l/2'S-i_i /o  ^   '^■>     "tl^en  using  the  definition 
(3.11a)  and  the  relation     min   (a,b)   =    j  [(a+b)  -   |a-b|] 


we  eet 


i'\^  '^  1-1 


gj    =     f   ^Sj-1/2   ^  gj+1/2   -   "j+l/2lgj+l/2   -  ^3-1/2' 

(3.13a) 

=  ^±1/2  ""    ^   ''*^gj+l/2  -  ^j-1/2^   -  ^j+l/2'^j+l/2  -  §3-1/2'^- 


Frcm  (3.11b)  we  conclude  that  if     v     is  Emooth  and     Q(x)     at  least  Lipschit: 
continuous  then 


'V  'V'  /    2  V 


Thus   (3.13a)  ar.d  (3.13b)   inply  that 

gj    =  gj,i/2  +  0(A^)-  (3.13c) 

It  is  easy  to  see  that  (3.13c)  holds  even  if  g--,_i/2  ■S4+1/2  ^  '^'  ^°^   ^^^^ 
g .  =  0  but  g-+2/9  =  O(A^)  itself  (since  A.  ^^^  =  0(A^)). 

Relations  (3.9a)  -  (3.9b)  follow  iinrrediately  by  rewriting  (3.13c)  as 


gj  =  ^,1/2  ^  0(a2)  ,       g.,,  ^  lj,,/2  ^  O^^'^' 


Ke  turn  nCTv  to  Dixsve  (3.12).  V'e  observe  from  the  definition  (3.11a)  that 
g.  and  g._|_,   cannot  be  of  different  sign,  hence 
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Igj+1    -    gjl     ^   "Bx(|gjl,|gj^^|) 

^     max[Min(Igj_^/2l'lgj+l/2l^'  "^^  lgj+1/2' ' '§3+3/2' ^^ 
-^     |gj+l/2l- 


Thus  it  follows  from  (3. lib)  that 


l^j+l/2l    =    Igj+1  -  gj'/|^j+l/2  ^'    ^^    'gj+l/2l/|    ^j+1/2  ^1 


I^S+1/2'    "   ^^j+1/2^^1      ' 


this  coirpletes  the  proof  of  C3.12). 

We  show  new  that  the  5-point  second  order  accurate  scheme   (3.8)  with  (3.11) 
is     TVNI     under'  the  same     CFL  restriction  of  the  original  3-point  first  order- 
accurate     TVNI      scheme   (3.1). 

LeKrra  3.^:      Suppose       Q(x)      satisfies   (3.2)  and  g.      is  defined  by  (3.11),   then 
the  finite-difference  scheme   (1.M-)  with  the  numerical  flux  (3.8)   is  TVNI     under 
the     CFL     restriction   (3.3). 

Proof:      Since    (3.8)   is   (3.1)  applied  to  a  modified  flux     f^      (3.8a),   it  can 

—  _M 

be  re'-NT-itten  as   (3.5)  withv.     ,  ,2    replaced   by         ^•+1/2      (3-8b).     We    conclude 

from   Lemma    3-1    that    the    scheme     (3.8)     is    TVNI    under    the   modified   CFL 

restriction 
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(3.14) 
^'3'"j+l/2l    "   " 

To  corrolete  the  nroof  of  Lemma  3.4  we  show  that   (3.14)   is    implied    by    the   original 
CFL     condition  (3.3).      UsLng  (3.12)  and   (3,2)  we  get 

l>l/2l    =    l^j-H/2  ^  ^j+l/2l    ^^    1^3+1/2'    '    1^3+1/21    ^    '^j+1/2'    + 
.1-1      (|v.,,/2l    -1)^<1 


vdienever 


|v-,-| /oj  5:  Q(v.  -,  ,^)  $  1;   this  shows  that  (3.14)  holds, 


Remarks:  (1)  If  A.  .,  ,„v  =  0  then  it  follows  immediately  from  the  definition  (3.11) 

that  £.  =  c.,.  -  0.  This  shcvjs  that  the  modified  numerici.1  flux  (3.S<i)  is  consistent 
3  D+1 

with  the  physical  flux  f(u)  Ln  the  sense  of  (1.4c).  Th.e  scheme  (3.8)  +  (3.11)  is  TVNI 
and  ■  therefore  it  has  convergent  subsequences  for  all  initial  data  of  bounded   .  ".  • 
^otal  va;r-iaticn:  the  limits  of  thsse  subsequences  are  we2"<  solutions  of  the  scalar 
conservation  law  (2.1).  To  coirplete  the  convergence  prco:  one  has  to  show  that 
all  these  limits  are  the  saire.   In  the  constant  coefficient  case  the  solution  to 
the  Pi'P  (2.1)  is  linique  arid  therefore  the  scheme  is  crnvergent  (note  that 
the  scheire  is  nonlinear  even  in  the  constant  coefficient  easel)  In  the  nonlinear 
case,  ccnvergence  ^sdll  follccvz  if  one  shows  that  the  scher?  is  consistent  with 
Oleinik's  entropy  condition  in  the  sense  of  (1.5).  We  srj^ll  discuss  consistency 
with  the  entix^py  condition  iii  Section  5 . 
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(2)  Condition  (3.6)  is  only  a  necessary  condition  for  second  order  accuracy. 
It  becanes  a  sufficient  condition  if  the  coefficient  in  the  0(A  )  term  in  (3.6) 
is  different iable ,  except  possibly  at  a  firjLte  nunnber  of  points  N(t,  At),  such 
that  NAt  ->  0  as  At  ->  0  for  all  t.   It  is  clear  from  (3.13)  and  (3.9)  that 

the  troublescine  points  where  the  scheme  (3.8)  +  (3.11)  may  degenerate  locally  to 

2 

OCa  )  truncation  error  are  those  vjhere  s.^,  ,„  (3.11c)  is  discontinuous,  i.e. 

2 
^ere  Q(v)  -v  =0   oru^=0.  The  fact  that  the  scheme  is  TV14I  controls 

the  possible  increase  of  the  number  of  local  extremum  points  in  the  conputed 

solution.  The  schen^.es  that  v;e  consider  in  section  5  all  have  the  monotonicity 

piroperty  (see  section  2);  i.e.  the  number  of  local  extremum  points  in  the  conputed 

solution  is  nonincreasing  in  time,  and  thus  bounded  by  that  of  the  initial  data. 

(3)  The  modified  equation  of  the  scheme  (3.1),  i.e.,  the  equation  which  it  aoproxi- 
Jiates  to  second  order  accuracy,  is  (2.7)  with 

e(u,X)  =[q(v)  -  v^  -^2,     V  =  Xa(u)  (3.15a) 

We  rewrite  the  modified  ecuation  as 


u^  +{f-^Q(v)  -  v^]  u  }  =0,  (3.15b) 

t      2x2  XX' 


and  observe  that  g-,-,  /^   (3.11b),  and  consequently  g.   (3.11a),  is  an 
approxirration  to  the  term 


g«^Q(v)  -  v^]  u^  =  I    [Q(v)  -  v^]  (Ax  u^).  (3.15c) 
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Our  irethod  to  convert  a  first  order  accurate  TVNI  scheme  into  a  second  order 
accurate  TVNI  scheme  is  based  on  the  following  heuristic  argument:  The  first  order 
scheme  C3.1)  approxiriHtes 

u^  +  \.f-jg\  =   0  (3.16a) 

to  second  order  accuracy.   Therefore,  applying  the  same  scheme  to 

u^  +  [f  +  ^]^  =  0  (3.16b) 

results  in  a  second  order  accurate  aporoxiination  to  u ,  +  f  =0.  To  ba  able 

t    X 

to  apply  the  scheme  to  (3.16b)  we  have  to  define  g(u)  with  a  bounded  derivative 
dg/du,  therefore  g-,-,  /^  (3.11b)  is  replaced  by  g.  (3.11a)  (see  [  s])- 
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U.   Syst&T5  of  Conservation  Laws. 

In  this  section  we  describe  how  to  extend  our  new  scalar  scheme  of  section  3 
to  systens  of  conservation  laws.  Our  extension  technique  is  a  somewhat  generalized 
version  of  the  procedure  suggested  by  P.  Roe  in  [  lU] .  The  basic  idea  is  to  extend 
the  scalar  scheme  to  the  system  case  by  apply iny  it" scalar ly"  to  each  of  the  (approp- 
riately linearized)  characteristic  variables. 

Let  S(u)  =  (R-'-(u),  R^(u),...,K^(u))  (4.1a) 

be  a  matrix,  the  coloms  of  vdiich  are  the  right  eigenvectors  of  the  Jacob ian 
natrix       A(u)    (1.2a).     Then 

o-Iao         .  -A  ^f    \x:  (U.lb) 

SAS  =  A  ,  A..=a(u)<5.. 

The  rows  L  Cu) ,  L  (u) , . . .  ,L™Cu)  of  S~  (u)  constitute  a  conplete  system  of 
left  eigenvectors  of  A(u)  which  is  bi-orthononral  to  the  system  of  right  eigen- 
vectors, i.e. 


L^R^  =6...  (4.1c) 

11 


In  the  constant  coefficient  case  A(u)  h  A  =  const. 


u^  +  Au  =  0,  uCx,0)  =  i|.(x),        -co  <  X  <  «•.         .      (4.2) 
t     X    '      ' 


k 
one  defines  characteristic  variables  w  =  (w  )  by 


w'^  =  L^'u,      w  =  S  ^u.  (4.3) 
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It  follows  from  (u.l)  that  (U.2)  decouples  into  m  scalar  characteristic 
equations,  1  ^  k  ^  m 

vt  +  a^  ^  =   0,  w'^(x,0)  =  L%Cx),       -co  <  X  <  ».         C+.U) 

t       X    ' 

This  offers  a  natural  way  of  extending  a  scalar  scheme  to  a  constant  coefficient 
system  of  equations  Cu.2)  by  applying  it  "scalar ly"  to  each  of  the  m  scalar 
characteristic  equations  C4.4). 


The  characteristic  variables  w   in  CU .  3 )  can  also  be  viewed  as  the  com- 
ponents  of  u  in  the  coordinate  system  {R  },  i.e. 


u  =  I   wV  (4.5) 

k^^l 

We  use  this  interpi-etaticn  of  chars cteristi^c  variables  to  extend  the  scalar 
scheme  to  general  nonlinear  systems  of  conservation  laws. 


Let  '^^+-^/2  -   V(v. ,  v._^^)  be  an  average  of  v.  and  v.^^,  i.e.  a  smooth 
function  V(u,v)  such  that 


V(u,v)  =  VCv,u)  (U.6a) 

V(u,u)  =  u  (4.6b) 
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and  let    a.^,  /2     denote  the  coiiponent  of    ^^+3^/2^  "  "^^+1  -  v.     in  the  coordinate 
system    ^^^^.^^^2^^ 

':n/2-  =  Jl''j*l/2  "^.1/2  "'•'^' 

°j+i/2  =  ''jn/2  'inn-"  -  <'*-™ 

here  v/e  use  the  notation  convention     ^A+1/2  ~  ■^^^i+l/2^  "  ^^^^"^'^  »^-i+i^^  • 

Vfe  now  extend  the  scalar  scha-ne  (3.8)   +   C3.ll)  to  general  systems  of 
conservation  laws  as  follows: 

f'  =  -?-«"f^*l/2-"f?-l/2''  <^-^-' 


fjn/2=    |C«"j'   +«-3+l>]  + 


(4.8b) 


^  ir   J,  ^j.1/2^^'  ^  ^>i  -  ^^>i/2  ^  Vi/2^  Vi/2^ 

where       v^^^^^  =  ^^^^^j+l/2^     ^^ 

4  =   s^^^^,2'^[0,  min(lSi+i/2Ugi_i/2-Si^l/2^^'      ^in/2  "   sign(g.^^/2^    (14.8c) 


W2   4      tQ^vJ,^/2>   -    ^Vl/2^'^  Vl/2  ^^-^^^ 

k        f^i+1  -  Si)/«i+i/2       ^^^^"     "in/2  ^  0 

^i+1/2.1  ^  (U.8e) 

(  0  vdien       ^^+1/2   ~   ^ 
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The  second  order  accurate  one- step  Lax-Wendroff  scheme  can  be  represented 


'j"i/2  =  I  [«-l>  *  «^m'1  -  k  J/"3n/2''  °^i/2  '^n/2   '^•^> 


.  1  :f(v  )  *   f(v  !>]  -  ^  [A(v.,^/,)j2s.,,/2V 


Lenma  4.1.   Siippose  {Q  Cx)}   are  Lipschitz  continuous,  then  (4.8)  satisfies 


f .   ,  =  f^"  ,  +  OCA^)  (4.10) 

3+1/2   ^j+1/2 


Proof:  Rev.T'ita  (4. 8b)  as 


where 


h   ^Q^^n/2  ^  >l/2^  -  Q^^Vl/2^^  Vl/2 


(4.11b) 


and  then  use  (3.10c)  and  conclude  (3.9)  from  Lemna  3.3. 

We  define  the  total  variation  TV(v)  of  the  vector  mesh  function  v  to 

be 

»    m   1 

j=-«'  k=l  -^ 

vtiere     a._|_    ,_      is  defined  by   (4.7),  and  shov;: 

Ler?ra  4.2     Suppose     Q^(x)     satisfies   (3.2)   for  all     k,     and  that     A(u)  =  A  =  constant 

then  the  scheire   (4.3)   is     Vjl^l     uiider  the     CFL     restriction 
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XnBxlA^I  $  V  -  Tidn  v^  i     1  (4.13) 


k 
whea^e  y   are  the  restrictions  in  (3.2) 


k      k         k 
Proof;   Because  of  the  assunption  A(u)  h  const,   {R  }  ,   {L  }  and  {a  } 

are  ^11  constant.  1-M.tiplying  (U.Sb)  from  the  left  by  L  ,  we  obtain  (3.8d) 

k         k       k 
for  the  characteristic  variable  w   (4.3);  g.  and  ^-^.-i/o  ^   (4.8)  becomes 

identical  with  (3.11).  Thus  by  Lemma  3.4  we  conclude  that  under  condition  (4.3) 

the  total  variation  of  each  of  the  characteristic  variables  is  nonincreasing, 

^nd  therefore  the  total  variation  (4,12)  is  nonincreasinp  as  well. 


Corollary  4.3.  The  scheme  (4.8)  in  the  constant  coefficient  case  is 
convergent  under  the  restriction  (4.13)  for  all  initial  data  of  bounded  total 
variation,  and  is  second  under  accurate. 

Vfe  remai-^k  that  this  corollary  is  not  trivial  since  the  schaTie  is  highly 
nonlinear  even  in  the  constant  coefficieiit  case. 

Our  technique  to  extend  scalar  schemes  to  the  system  case  does  not  require 
any  particualr  form,  of  averaging  V(u,v)  (4.6)-  Roe  in  P-5  ]  is  using  a  specific 
fonn  of  averaging  that  on  top  of  being  mathenHtically  pleasing,  also  enables  the 
conputational  advantage  of  perfectly  resolving  stationary  discontinuities. 

In  [  7  ]  we  sha-;  that  if  the  system  of  conservation  laws  (1.1)  possesses 
an  entropy  function  (1.3),  then  it  is  symmetrizable ,  aj-id  there  exists  a  mean  value 
Jacob ian  A(u,v)  such  that 
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(i)  f(v)  -  f(u)  =  A(u,v)  (v-u)  (4.1Ua) 

(ii)  A(u,u)  =  A(u)  (U.lUb) 


(iii)  A(u,v)  has  real  eigenvalues  {a  '(u,v)}^_-|^  and  a  conplete  set  of  right 


ignevectors  {K  (u,v)},  _, 


eign' 


In  the  conte>rt  cf  the  scherre  (4.8)  Roe's  extension  technique  is  expressed  by 

k  k 

takin?  A. j,  ,„  and  R.,-  ,o  in  (U.7)  -  (4.8)  to  be  the  eignevalues  and  the  right 
j+1/2       j+i/2 

eigenvectors  of  the  mean  value  Jacobian  A(v.,  v._^,)  (U.14a),  respectively. 

]^ 
Thus  if  a-^-,  /^  are  defined  by  (4.7a) 


k       r-X 


(4.15a) 


V.,,  -  v.  =   >   Ct.,,  ,_  r<-,-.,r, 
^^    ^   k=l   ^  ^   '^    "-^^ 


then  it  f oilers  from  (4.14a)  that 

^   V      V      V 

fCv,,^)  -  fCv.)  =  J^  „^,,,2  aj,,/,  F^,,,j.  C.lSb) 

The  relation  (4.15)  makes  the  scherre  (4.8)  a  more  faithfu!!.  extension  of  (3.8) 
in  the  sense  that  (4.8)  for  Tn=l  is  identical  with  the  scalar  scheme  (3.8). 

Wa  observe  th^t  if   f(u)  =  f(v)  in  (4.14a)  then  v  -  u  is  a  right  eigenvector 

or  A(u,v) ,  cc'-'^socnding  to  a  zero  eignvalue.  Hence  in  (4.15)  a  (u,v)  =  0  for 

k  i^  k   ajid  a  °(u,v)  =  0,  for  some  k  .  It  is  easy  to  see  that  if  Q  (0)  =  0 
o  o 

in  (4.8)  Then  the  stationary  discontinuity 

u     x  <  0 
u(x,t)  i  d)  (x)  =  {  ,  f  (u)  =  f  (v) 

V     x  >  0 
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is  also  a  stationary  solution  of  (U.8).  (see  [15]). 

In  the  case  of  the  Euler  equations  of  gasdynamics,  where  the  flux  f(u) 
is  an  hcxiTogeneous  function  of  u  of  degree  1,  it  is  possible  to  express  A(u,v) 
in  C4.1U)  as 

ACu,v)  =  A(VCu,v)).  (4.17) 

This  relatively  sinple  function  V(u,v)   (see  [15])  will  be  described  in 
Section  7. 

Renarks:   (1)  Note  that  we  use  0  (x) ,  thus  allov;ir.g  different  functions  (3.2) 
for  different  characteristic  fields.  As  observed  by  P.  Roe  [14j  the  extension 
technique  of  this  section  permits  even  the  use  of  conpletely  different  scalar 
schemes  for  different  characteristic  fields. 

(2)  In  most  applications  it  iiHy  be  advantagous  to  replace  the  term 

i  [f(v.)  +  f(v.^,  )]  in   CU.Bb)  by  fiv.^^,^)   =   f(V(v.,  v.^, )).  For  sur^  this 
2      3       :+l  ^         j+l/z  3'  3+1 

simplifies  the  progrBinniing  and  reduces  the  CPU  time,  without  altering  the  main 
propearties  of  the  scheme.  A  possible  disadvantage  of  such  a  change  is  that  the 
scheme  nay  lose  the  property  of  perfect  resolution  of  stationary  discontinuities. 

(3)  The  particular  definition  (4.12)  of  total  variation  is  motivated  by  the 
definition  of  Glirrm's  functional  in  [3  ].  When  applied  to  a  piecewise- smooth 
solution  u(x,t)  of  (1.1) 

"   m   ■,  m    , 

lim  ■IV(u)  =  /   n^-  (u>Uvl^'^  ^  ^   ^   l«^x.)|  (4.18) 

AX-K3       -o"  1<f1  i  ]<f1      ^ 
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k  k 

where  x.  are  points  of  discontinuity,  and  a  (x. )  denotes  the  value  of  "^j^+t/o 

in  C4.7)  evaluated  v;ith  respect  to  v .  =  u  ((x.  )_,t),  v^^.^^  =  uC(x.  )_|_,t) . 


There  is  no  reason  to  expect  that  the  functional  (U.18),  and  consequently 
C4.12),  is  generally  nonincreasing  with  t.   Based  on  Glimm's  results  [3  ]  we 
do  hov;ever  believe  that  this  functional  (under  certain  conditions)  is  bounded  in 
t.  At  this  tine  v;e  do  not  hs-ve   estirrates  of  the  possible  increase  in  total 
variation  in  solutions  of  the  scheme  (U.8),  and  therefore  cannot  prove  convergence 
in  the  nonlinear  system  case. 
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5.   On  the  choice  of  0(x). 

In  section  3  we  have  presented  the  basic  idea  of  our  new  schejne  in  the 
scalar  case,  v;hich  can  be  described  algorithmically  by:  "Take  a  3-point  TVNI 
schene  (3.1),  which  is  first  order  accurate,  and  apply  it  to  the  modified  flux 
value  (3.3)  to  obtain  a  5-point  second  order  accurate  TVNI  scheme  (3.11)  +  (3.8) 

Vfe  begin  section  5  by  considering  different  choices  of  QCx)  in  the 
3-point  scheme  (3.1).  In  the  end  of  this  section  we  argue  that  most  of  the  pro- 
perties of  the  3-point  point  scheme  (3.1)  go  over  to  the  mxiified  5-foirit  scherK 
G.8)  +(3.11). 

In  section  6  we  discuss  the  system  case. 

A  natural  cchice  of  Q(x)  .in  (3.2)  is  Q(x)  =  |x|  ,  which  corresponds  to 
the  least  dissipative  TVl^II  scheme  of  the  form  (3.1),  Itie  scheme  (3.5)  with 
Q(x)  =  jx|  can  be  re^^fritten  as 


n+1    n   ,-     V-  .      n       f-  v+  .      n  /-  n  n 


v^ere 


V  =  min(v,0)  =  y  (v  -  |v|),  v   =  nHx(v,0)  =  y  (v  +  Iv|).       (5.1b) 


This  sche^fi  (5.1)  is  a  generalization  of  the  well-known  upstream 
differencing  sdneme  of  Courant,  Isaacson  and  Rees,  and  it  is  well  investigated 
in  the  literat'ure  (see  [5],   [10]  and  the  references  cited  there).  We  now  give 
a  brief  review  of  its  relevant  properties. 
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Leirrra  5.1.      Solutions  of  the  upstream  differencing  scheire   (5.1),  under  the     CFL 
restriction 


Xrrax   |a(u)|    <   1,  <i>  =  range  of  initial  data,  '  (5.2) 

ueo 

have  the  nono tonicity  property    stated    at    the    beginning   of    Section    2. 


Sketch  of  Proof:      Let     \Kx)     be  the  piecewise-linear  interpolant  of     v   ,     i.e. 

V/(x.)   =  V.    ,     X.    =   jAx. 
3  3  D        ■" 


Then  (5.1),   under  the     CFL     restriction  (5.2),   ijrolies  that  there  exists 
a  Tnonotone  ncndecreasing  sequence  {  x . }      such  that 

v^^     =  VKx.),  I  X.   -  X. I    ^  Ax 

^  ^  3  3 

This  sho".-7S  that  v  '    has  no  rore  local  extremum  points  then  v  ,  and  that 
the  value  of  a  local  minirrum  (raxai'TijjTi.)  is  nondecreasLng  (nonincreasing) ,  which 
is  the  assertion  of  Lenrva  5.1. 


We  dbsejrve   that  the  stationary  iunp  discontinuity  (U.16),  admissible  or 
inadmissible,  is  also  a  satationary  solution  of  the  upstream  differecing  scheme 
(5.1).   Cn  one  hand  this  is  a  desirable  property  as  it  inrplies  good  resolution  of 
staticrer^'.'  shocio ;  an  the  other  hand  it  indicates  that  the  scheme  many  select 
ror:r:hvsical  v:ea!<  solutions  that  do  not  satisfy  the  entropy  condition  (1.3b).  This 
property  is  related  to  the  fact  that  the  viscosity  term  B(u,X)   (3.15a)  vanishes 
for  V  =  Q,  and  it  is  ccs?r?or\   to  all  schemes  (3.1)  with  Q(C)  =  0.' 
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To  prevent  entropy  condition  violation  of  this  sort  and  wake     Q(x) 
smoother  at  the  same  tims,  we  modify  Q(x)  =  |x|  near  x  =  0  to  be 


QCx) 


^^e  for   |xU2e  ^^^3^ 

.1x1  for   lx|  5  2e 


with  say  e  =  0.05.   Csee  Fig.  1) 

This  change  has  increased  the  amount  of  nurrerical  viscosity'  for  |x|  <  2g 
so  that  now  6(u,A)  >  0  for  \\>\    <   1.   6(u,X)  vanishes  only  for  v  =  1; 
which  can  be  handled  by  taJcing  u  <  1,  say  p  =  0.95,  in  the  CFL  restriction 
(3.3).  Consequently  the  scheme  (3.1)  with  Q(x)  (5.3),  cannot  have  perfect 
'resolution  of  the  stationary  shock  (4.16). 


Another  possible  choice  of  Q(x)  which  correspond?  to  a  schei.-ie  that  has 
aljpeady  been  investigated  is  Q(x)=x  +-j-  (see  L  ^-1     -".nd  [5]).  This  scheirie 

is  the  lax-Wendroff  scheme  modified  by  the  addition  of  "::'.'.3  linear  viscosity  term 

1  12 

■5<v._|_-  -  2v.  +  v.  ^).  Since  now  QCx)  -  ]x|  =  ([x]  -  -)  ^0,  it  is  the  require- 
ment Q(x)<:.l  in  (3.2)  which  restricts  the  CFL  condition  (3.3)  to 
V   =  ^  /2»  0.866  <  1  (see  Fig.  1). 


We  turn  new  to  discuss  the  second  order  accurate  TVNI  scheme  (3.8)  +  (3.11). 
We  note  that  the  truncation  error  of  this  scheme  is  OCJ.  ).  except  possibley  where 
Ujj  =  o,  independent  of  the  particular  form  of  Q(x).  Th.L:  the  modified  scheme 
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reJTDves  all  0(A)  errors  due  to  nuinerical  viscosity,  except  at  the  points  where 
u   goes  through  zero;  there  some  viscosity  is  needed  to  ensure  monotonicity . 
Consequently  we  expect  the  new  scheme  to  have  high  resolution  of  shocks,  stationary 
or  noving,  alirost  independently  of  the  particular  choice  of  Q(x). 

Eased  on  these  considerations  it  seems  to  us  that  C5.3)  is  the  best  choice 
of  QCx)  for  C3.ll) 

Ifext  we   deal  with  the  question  v.-hether  properties  of  the  schem.e  (3.1) 
go  over  to  the  irodified  scheme  (3.8)  +  (3.11).  We  expect  properties  of  the  origijial 
schejre  that  can  be  expressed  as  a  relation  bet'^een    F(v  )  and  F(v  '  ),  where 
F  is  soiTE  functional  that  depends  on  mesh  values  only,  to  go  over  to  the  modified 
second  order  accurate  scheme.  We  have  seen  th^'.t  this  is  true  for  the  TVTC  pro- 
perty; v.'e  can  also  prove  that  for  Q(x)  =  jx]   the  modified  scheme  has  the  mono- 
tonicity property,  and  th-.--t  in  thie  stationary  discontinuity  case  (4.18)  it  behaves 
the  sane  way  as  the  original  scheme. 

To  support  our  conjecture  in  general,  we  note  thaT  the  modified  scheme  is  the 
original  scheme  applied  to  a  rcdified  flux  f  +  j  g     (3.15b).  Thus  if  g  were 
strictly  a  function  of  u,  our  statement  becomes  trivially  true.  Hov;ever  g  in 
C3.ll),  which  is  consistent  with  zero  flux,  rray  be  cosidered  g(v  )  only  where 
v^  is  ironotone.  The  particular  form  used  in  (3.11)  is  designed  to  match  the 
different  definitions  of  gCv")  in  neighbo'jring  monotone  sections  by  setting 
g  =  0  at  local  extremum  poLnts  (see  [  5]).  That  this  matching  is  smooth,  is 
evident  by  the  fact  that  y  -   Ag/iu  is  well-defined  eveiY.-;here ,  and  is  bounded  in 
absolute  value  by  (3.12). 
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Hopefully  the  piecewise  definition  of  g(v  )  does  rrat  spoil  validity  of 
such  properties  for  the  modified  scheme. 

Based  on  these  arguments  v;e  form  the  conjecture  that  if  the  original 
scheme  C3.1)  is  consistent  with  the  entropy  inequality  (1.3b),  so  is  the  modified 
second  order  accurate  scheme  (3.8)  +  C3.ll). 
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6.    On  resolution  and  entropy. 

The  extension  technique  of  section  4  consists  of  applying  scalar  schemes, 
not  necessarily  the  sare,  to  characteristic  fields.  Therefore  we  1-iave  the  possibi- 
lity to  ousted  fit  the  scheme  to  the  conputational  needs  of  each  characteristic 
field, 

Vfe  consider  here  systems  of  conservation  laws  where  the  characteristic  fields 

are   either  gen-'oinely  norilinear  (a'^R  ^  0)     or  linearly  degenerate  (a  K^  =  0, 
— -■ u  ' ^ u 

see  [11]).  The  waves  of  a  genuinely  nonlinear  field  are  either  shocks  or  rarefac- 
tion waves,  depending  whether  the  characteristics  are  com/ergent  or  divergent. 
The  waves  of  a  linearly  degenerate  field  are  exclusively  ccntact  discontinuities. 

First  let:  us  consider  the  latter'  case  of  a  linearly  degenerate  characteristic 
field.   It  is  well  known  that  the  RiemarJi  invariants  of  tYls   field,  one  of  which 
is  the  characteristic  speed,  ai^  continuous  across  a  contact  discontinuity.  Tnere- 
fore  the  propagation  of  a  contact  discontinuity  is  ccnputa-ionally  equivalent 
to  that  of  a  scalar  discontinuity  moving  with  a  constant  characteristic  speed. 
As  remarked  in  section  3,  the  solution  of  xhe  first  order  accurate  scheme  (3.1) 
is  the  saire  ud  to  second  order  terms,  as  the  solution  of  the  modified  differen- 
tial equation  <3.15b)  Cor  C3-15a)).  The  modification  of  the  flux  f(u)  =  constant. u 
by  the  addition  of  the  terra  C-  ^  g)  has  the  effect  of  nvijcing  the  cnaracteristic 
field  slightly  divergent  in  a  region  of  a  discontinuity.  Consequently  the  computed 
discontin^jity  is  being  spread  at  the  rate  of  y^,     where  n  is  the  number  of 
time-steps  taken;  similarly,  a  p-th  order  accurate  standard  scheme  spreads  a  contact 
discontinuity  a1:  rhe  rate  of  n   ^     (see  [5  J). 
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Our  new  second  order  accurate  scheme  is  obtained  by  applying  the  original  one 
to  the  flux  f  +  J  g     (3.16b).  This  has  the  effect  of  reducing  the  divergence  of 
the  characteristics  induced  by  the  original  scheine  in  a  region  of  discontinuity, 
and  we  expect  the  corputed  discount inuity  to  spread  li>;e  n   ,  rather  than  n 


To  ensure  that  the  conputed  contact  discontinuity  does  not.  spread  in  time  at 
all,  we  want  to  make  the  effective  characteristic  field  slightly  convergent.  To 
acccnplish  that  we  increase  the  size  of  g  in  (3.16b)  e.g.  by  multiplying  the  RHS 
of  C4.8cl  by  a  +  e^),  i.e. 

g^  =  (  1  +  ei)sJ^^/2  "^l^^'  "^   ^l%+i/2''Si-l/2-^i+l/2^^         ^^-^^^ 
where 

a._^,  ,^  is  defined  by  (4.7). 


In  £6  J  we  show  that 


C  $  6^  ^  1  (6.2a) 


e.  =  OCax)  in  regions  of  smoothness.  (6.2b) 


and  that  6-  =  0(1)  in  regions  of  discontinuities. 


It  is  easy  to  see  that  this  change  in  (4.8c)  modifies  the  numerical  flux 

2 

(4.8b)  by  0(a  )  and  thus  does  not  spoil  the  second  order  accuracy  of  the  schene. 

Frxm  (6.2a)  and(3.12)  it  follo\vrs  that  now 
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l^?.l/2l  -  ISj-l  -  4'^'Vl/2'  '\^^^U/2^  -   ^^jn/2^  I         ^'-'-^ 
Hence  for  Q(x)  (3.2)  we  get  in  Lenma  3.4 


I^j+l/2l  "^Iyj+1/21  ^  l^j+1/2'  ""  Q^^j+1/2^  -  ^^j+1/2^  ^^'^^^ 


For  QCx)  =  |x|  we  have  that  the  RHS  of  (6.3b)  is  less  or  equal  1  for 
\^'+'\/?\    ^   ■'"'  ^^"^  consequently  the  modification  (6.1)  does  not  alter  the  CFL 

condition  for  this  field.  The  same  is  true  for  Q(x)  in  (5.3)  with  e<  1/3. 

2 
(Hcfi-jever  for  Q(x)  =  x  +  1/4  (6.3b)  implies  a  more  stringent  restriction, 

naively     max|v^_^  .^l  ^  3/4). 

We  reirark  that  there  are  no  entropy  considerations  associated  with  a 
linearly  degenerate     characteristic  field. 


We  turn  now  to  the  case  of  a  genuinely  nonlinear  characteristic  field, 
where  the  corputational  aspects  to  be  considered  are  resolution  of  shocks  and 
enforcei!>ent  of  the  entropy  condition. 

Unlike    contact   discontinuities,    shocks    are    formed    and   main- 
tained   by    local    convergence    of    the    characteristic    curves.       The 
reduction    in    numerical    viscosity    due    to    the    addition    of     (4.8c) 
is    usually    sufficient    to    ensure    good   resolution  of    shocks.       When 
higher    resolution    is    desired,    one   may    employ   a   mechanism  of    the 
form    (6.1)     to    enhance    the    local    convergence    of   the    characteristics 
in    the    shock    region    and    thus    improve    its    resolution    (see     [5]). 
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The  conputational  aspect  of  the  entropy  condition  nay  be  demonstrated  by 
considering  a  RieiiHnn  problem  for  (1.1)  involving  two  vearhy  states  u,  and  Un. 
Lax  in  EllJ  proves  (assuming  all  fields  to  be  genuinely  nonlinear)  that  there 
exists  a  unique  "physical"  solution  to  this  problem,  vMch  consists  in  general  ■ 
of  both  shocks  and  rarefaction  waves.  One  can  show,  using  exactly  the  same  method 
of  proof,  that  there  also  exists  a  weak  solution  consisting  only  of  shocks  ,  which 
may  differ  from  the  physical  one.  Whenever  the  first  contains  rarefaction  waves, 
the  letter  ]-ias  nonphysical  "expansion  shocks".  Therefore  to  ensure  that  the  scheme 
selects  a  physically  relevant  weak  solution,  consistency  with  the  entropy  condition, 
on  top  of  consistency  with  (1.1),  is  required. 

In  section  5  we  have  argued  that  the  scalar  scheira  (3.8)  +  (3.11)  is  consistent 
with  Oleinik's  entropy  condition.  The  question  arises  now  whether  the  system 
version  (U.8)  of  this  scheme  will  always  select  the  ph^^ical  weak  solution, 
T.P,Liu  ll2j  shCT^7s  that  if  a  weak  solution  of  (1.1)  con-£U.ns  only  "admissible" 
discontinuities  that  satisfy  a  "scalar"  Oleinik  conditio.  v;ith  respect  to  the 
Rankine-Hungcniot  curve,  then  this  solution  also  satisfi-=  the  entropy  inequality 
(1.3b).  Using  this  theory  the  above  question  can  be  reformulated  as:  Will  the 
Scheme  C^^.S)  reject  discontin'uities  that  do  not  satisfy  Tleinik's  condition  along 
the  Rankine-Hugoriot  curve?  It  seems  possible  to  answer  this  question  affinratively 
for  sufficiently  weak  shocks ,  and  feasibly  for  a  larger  rlass  of  problems . 


1  We  reiiark  that  there  is  also  a  one-parameter  family  zf  states  connecting  u, 
and  Uq  through  rarefaction  curves  only  (see  D-3  J).  R.v.ever  a  "negative  rare- 
faction U'ave"  is  not  carputationally  realizable  as  it  is  a  multi-valued  function 
of  X. 
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7.    Application  to  Euler  equations  of  gasdynandcs 

In  this  section  we  describe  the  application  of  our  new  scheme  (4.8)  to 
the  Euler  equation  of  gasdynamics: 


w^  +  f(w)  =  0 


(7.1a) 


(P    ) 


f(w)  =  uw  + 


^  1 


(7.1b) 


p  =  (y  -  1)  (E  -  1/2PU  ), 


(7.1c) 


Piere     p,  u,  p     and     E     are  the  density,  velocity,  pressure  and  total  energy, 
respectively;     m  =     pu     is  the  inomentum  and  we  take     y  =  1.4. 


The  eigenvalues  of  the  Jacobian  iiHtrix    ACw)   =  f      are 


a,Cw)  =  u  -  c  ,  a2(w)  =  u  ,  a.^(.vi)   =  u  +  c 


(7.2a) 


1/9 

where  c  is  the  sound  speed,  c  =  C  yp/ p) 


The  corresponding  right  eigenvectors  are 


K  .  Cw^)  =    u  -  c 
1 

H  -  uc 


f   1) 


R,(w) 


1  2 
±   u 
2 


I.  o  ^  J 


,   R3 (w) 


1 

u  +  c 
H  +  uc 


(7.2b) 
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where  H  =  (E  +  p)/p  =  c^/Cy-l)  +  7  ^^  is  the  enthalpy. 


Let  a  (.w,  ,vi^ )  ,  k  =  1 , 2 , 3 ,  be  the  solution  of  the  following  system 
of  linear  eauations  (4.7) 


3     , 
w„  -  w,  =  I     a^^CVCw,  ,w^)).  (7.3a) 

l\       Li     1   1  Li    K 

k=l 

where  V(wj^,  w_)  (4.6)  is  some  average  state;  denote  its  velocity  and  sound 

-      -  k 

speed  by   u  and  c,  respectively.  To  calculate  a   in  (7.3a)  we  first 

evaluate 


(\   =  (y-D  {[E]  +  I   u^Ep]  -  u  [m]}/c^  (7.3b) 


C2  =  {[m]  -  urp]}/£  (7.3c) 


where  [b3  denotes  Cb]  =  b„  -  b,  ;  then  a   in  (7.3a)  are  obtained  by 


i""-  =  ^  (C^  -  C^),   a^  =  [p]  -  C^,   a^  =  ^  (C,  +  C^).  '        (7.3d) 


The  second  characteristic  field  corresponding  to  the  eigenvalue  u  is 

2     2 

U-nearly  degenerate,      i.e.,     a     R     h  0:     The  other  ch?racteristic  fields  corres- 

w 

ponding  to  the  eigenvalues       u  ±  c     are  genuinely  nonlliear.     Therefore  in 


2  13 

(6. la)     we  define     e._^,  ,„     by   (6.1b)     for     k  =   2     and  set     Q-^./o   =   9-_. 
j+1/2       -'  j+1/2         3  + 


1/2  =■  °- 


Given  wv  we  now  list  the  ooerations  needed  to  calculate  w.   : 

3  '  : 
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(i)       Ccopute:    (u-    =  Ta./P-)     and     (p.)     by     (7.1c) 


(ii)     Coirpute     u-^-./^     and     c.^-,  ,«     from     V(v.,v.^-,);     calculate 


ITHX. 

3 


^1^3+1/2'    ""  ^j+1/2^'      evaluate     a?^^^^'      ^'  =   ^'^'^     by   (7.3b)   -   7.3d), 


Define     X    =  At/Ax  =  v/max.  (I  u.^,  ,„|    +  c.  ,-,  z^)     where     u      is  the  prescribed 

j  '  j+1/2'    i+1/2  '■ 

CFL  restriction  in  (4.13) 


(iii)   Compute:  ^^^^^2   '  ^^"j  +  1/2  ^^  (7.2a);   g.^^/^  ^^'  ^^-^d); 

^j+1/2  ^^^  ^^'-^^^  ^°^  ^  ^  2'   "'^^  'j+1/2  ^  ^^n/2  =  °- 


(iv)  Calc-olate  g^  by  (6.1a). 

(v)   Compute  Y^+2/2  ^^  (4-8e),  and  "•+-^/o  by  (4.8b)  and  relations  (7.2b) 


(vi)  Corrpute  v/'"*''^  by  (4.8a). 


We  note  that  in  ccrnnutijig     f.  .^  ,„     in   (v),   one  could  take  advantaf!;e  of  the 

]+x/z 

k 
sinple  form  of  the     R       in  (7.2b). 


Nex"t  we  show  how  to  iiiplemeiit     Roe's     linearization  technique      (4.14)   - 
C4.J5)     ixi  the  above  algorithm.      Roe     presents  a  pajr^ticular  form  of  averaging 
VCw.  ,w„)      such  that   for  the  Euler  equations  of  gasdyric'LTacs,   the  iiean  value 
Jacobian  A(w,  ,w„)      :n   (4.14a)   can  be  e;<pressed  bv   (^.17).     Thds  averaging  talcos 
the  follcwii-ig  foiTn: 
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yl/2    N//l/2>  A  /1/2„>,^1/2S 

(7.Ua) 


^j+1/2  =    ^^^-^^    ^^j+1/2  -       l^j+1/2^^^^^ 


where  <'b>  denotes  the  arithemetic  mean 


<^b>  =  i  (b.  +  b.^,)  (7.4b) 


Therefore  to  use  Roe's  linearization  in  our  scheme  all  one  has  to  do  is  to 
coipute  u.^^^2  arid  c._^^^^     in  (ii)  by  (7.4). 

Vte  reiTork  that  the  averaging  in  (J .^)   is  rather  expensive.  It  seems  to  us 
that  in  nany  applications  the  sinple  arithmethic  average  (7.4b)  will  do  just  as 
well. 
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8.  ?h-iiT-firical  experiments. 

In  this  section  vre  present  some  numerical  experiments  that  demonstrate 
the  perforrrance  of  the  proposed  second  order  accurate  scheme.     We     consider  here 
the  follcwing  two  versions  of  it: 

1  1      ^        V  V 

^jn/2  =      2   tf(Vj)   +  f(v.,,)  -     ^- J^  6j,,/2  ^>l/2^  (8-^) 


where 

^j+1/2   =   ^   ^^^j+1/2   "  ^j+l/2^"j+l/2   -    ^Sj   ^  gj+i 


-   0^^4+1/2   +  y'^+1/2 ^4+1/2   -    ^4  +  4+1 )  (8-2^ 


and     g.      is  defined  by    an    appropriate    variant    of     (6.1)    and     (4.8d) 

)^ 
Y-,-.  /p    is    defined    by     (4.8e), 

The  first  version  is   (8.1)  -    (8.2)  \-n.th 


Q^(x)   H    |xl  (8.3a) 


and  u-lll  be  referred  to  as  the  scheme  ULTl;  The  second  version  is  (8.1)  - 
C8,2)  '^dth 

q'^(x)  =  x^   +   1/U  (8.3b) 

and  will  be  referred  to  as  the  scheme  ULT2. 

For  ccnparison  sa!<e  we  also  pi^sent  calculations  with  the  following  tv.-'o 
schaT.es : 
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(i)   The  second  order  accurate  Lax-Wendroff-type  scheme  (8.1)  with 


,k     _  ,  k    .2 


(8.4) 


which  is  referred  to  as  the  LW  scheme. 

(ii)  The  first  order  accurate  Godunov-type  scheme  of  Roe  (sne  [15]  ar.d  [10]) 
\^iich  is  defined  by  C8.1)  with 


^j+1/2  -   l''j+l/2' 


C8.5) 


and  is  referred  to  as  the  ROE  scheme. 

In  all  the  schemes  and  experiments  reported  herein  we  use  the  Roe 
linearization  C7.3)  -  (7.4). 

(I)   The  shock  tube  problem. 

Vfe  cosider  now  a  Rieirann  problem 


w(x,0)=. 


X  <  0 
X  >  0 


(8.6a) 


for  the  Euler  equations  of  a  polytropic  gas  C7.1).  Our  first  set  of  data  is 


\-- 


'0,U45 
0.3111 
[8.928  J 


0 

5 

0 

1 

4275 

(8.6b) 
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Other  numerical  experiments  vd.th  this  problem  are  reported  in  [6  ]  and  the 
references  cited  there. 

In  figures  2a,  2b  and  2c  we  show  the  results  obtained  by  the  ROE,  LW  and 
ULTl  schemes,  respectively.  The  numerical  values  are  shown  by  circles;  the 
exact  solution  is  shown  by  the  solid  line.  The  calculations  in  figure  2  were 
perforrned  with  100  time-steps  under  the  CFL  restriction  p  =  0.95  in  (4.13), 
and  140  cells. 

In  figure  3  we  repeat  the  calculations  presented  in  figure  2  for  a  diffei-ent 
set  of  data  for  the  Rieirann  Problem  (8.6a) 


\  - 


0 

125' 

0 

_  0 

25 

Other  numerical  experiments  ^^d.th  thi.s  problem  are  preserved  in  [16],  The 
calculations  in  figure  3  were  perforrried  with  50  tirr.e-steps  under  the  CFL  restriction 
P  =  0.95  in  (4.13),  with  100  cells. 

Ife  remark  that  the  solution  of  ULT2  for  these  probelms  look  almost  identical 
to  those  of  ULTl.   In  both  schemes  we  find  almost  no  dependence  on  the  CFL 
number. 


dl)  The  Q^jasi  1-D  nozzle  problem. 

We  consider  an  axisymnietric  nozzle  with  a  cross-section  area  A(x) .  The 
cross-secticn  average  of  the  flow  satisfy  the  following  one-dimensional  system 
of  equations 


w^  +  f(w)^  = 


-  1+6  - 


-    S(V7,X)  ,  S(W,X)     = 


■  0 

' 

p 

dA 
dx 

0 

- 

v^Tere     w,   f(w)     and    p     are  given  in  (7.1). 

2 

In  figures  4,5  and  6  we  present  numerical  approximations  to  steady  state 

solutions  of  (8.7a). 

In  figures  U  and  5  we  show  solutions  for  a  diver^gent  nozzle  vjith  the  cross- 
section  area. 


ACx)  =  1.398  +  0.3U7  tanh  (0.8x  -  4); 


(8.7b) 


the  flow  condition  is  supersonic  at  the  entrance  and  subsonic  at  the  exit. 
Figures  Ua  and  Ub  sho^-7  steady  state  solutions  on  a  crude  ;nesh  of  the  ROE  and 
ULTl  schemes,  respectively.  Figure  5  shows  the  ULTl  res'jlts  for  the  same  problem 
on  a  finer  mesh. 

In  figure  6  we  shc^-i   a  steady  state  solution  of  z'r.e.   'XTl  scheme  for  a 
convergent -divergent  nozzle  with  the  cross-section  araa 


1  +  (Aq  -  1  )  (1  -  x/5)' 


A(x)  = 


X  <^  5 


1  +  (A^  -  1)  [(x  -  5)/(x^  -  5)]^     X  > 


(8.7c) 


where  A  =  entrance  area,  A_,  =  exit  area;  Here  the  flow  is  S'ibsonic  at  the 
entrance  as  well  as  at  the  exit. 


2    These  figures  ar^   by  courtesy  of  Helen  C.  Yee  of  the  KASA-Ames  Research  Center. 
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The  exact  solutiors  in  figures  4  to  6  are  shown  by  the  solid  curve;  the 
values  of  the  numerical  solutions  are  indicated  by  a  rcsnbus. 

Hopefully  a  detailed  report  of  these  calculations  and  the  particular 
appxxsach  to  steady  state  will  be  published  elsewhere. 

(Ill)  2-D  Flow  through  a  duct. 

In  figures  7a  and  7b  we  show  solutions  to  the  problem  of  the  flow  of  air 
through  a  duct  containing  a  step.   Initially  the  flow  is  everywhere  to  the  right 
at  Mach  3,  with  p  =  1.4,  p  =  1  and  c  =  1,  The  duct  width  is  1,  its  length  is  3, 
and  the  step  of  bdeght  0.2  is  located  a  distance  of  0.6  from  the  entrance. 
Figure  7  shows  the  results  at  t  =  U  with  a  crude  uniJoiTti  Cartesian  grid  with 
AX  =Ay  =  0.1. 

This  probla-n  was  used  by  V/oodward  and  Colella  to  test  the  performance  of  various 
nuiiHrical  scheirtes  (see  [18]  and  the  references  cited  hitere). 

The  solutions  in  figure  7  v;ere  obtained  by  a  Strang  -type  dimensional 
splitting  of  the  form 

v"^2  =   Lv"  (8.8a) 

L-L^L^L^l^  (8.8b) 


\stiere     L       and     L       are  one-dimenional  finite  difference  ocerators 
X  y 

approxir.k?.ting 
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Lj^  :  w^  +  f(w)^  =  0   ,    L  :  w^  +  g(w)  =  0.  (8.8c) 

If  L   and  L   are  stable  and  dissipative  second  order  accurate  approximations 
to  the  one-dimensional  equations  in  (8.8c),  then  the  schene  (8.8a)  -  (8.8b)  is 
a  stable  second  order  accurate  approximation  to  the  2-di3i!iensional  problem. 

w^  +  f(w)^  +  g(w)  =0.  (8.8d) 

In  figure  7a  we  show  for  ccsnparison  sake  the  results  of  the  second  order 
accurate  hybrid  schene  (see  [  6])  (8.1)  with 


"W  =  ^'^Un^'  *  ^  "Un^  "Im  "•^' 


where  ^^^^./l     ^^   (6. lb)  for  all  k. 


In  figure  7b  we  show  the  results  of  ULT2. 

Both  figures  7a  and  7b  shew  30  equally  spaced  density  contoijrs;  Both 
calCTilations  were  performed  with  a  CFL  restriction  of  C.75 
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Vfe  reiTBrk  that  the  comer  of  the  step  in  our  calculations  is  treated  as 
a  sharp  comer  without  any  rounding  (or  equivalent  addition  of  nuirerical  viscosity) . 
The  sonic  line  emerging  from  this  comer  is  a  curve  on  which  an  eigenvalue  of  (8.8d) 
vanishes.   It  is  interesting  to  note  that  the  results  of  the  Godunov  scheme  in 
figure  la  of  [18]  indicate  that  part  of  the  sonic  line  rray  turn  into  an  expansion 
shock.  Vfe  find  a  similar  behaviour  in  the  results  of  ULTl  for  this  problem. 

Altogether  V7e  find  the  performance  of  the  new  second  order  accurate  scheme  to 
be  quite  pleasirig.  V/e  note  that  the  scheme  is  simple  to  program  and  requires  only 
slightly  more  CPU  time  than  a  Lax-Wendroff  scheme  with  some  artificial  viscosity. 
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Figure  2a.    ROE  scherrie  for  (3.6b) 
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FigTjre  2b.    LW     schers  for  (S.6b) 
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Figure  2c.  ULTl  schenie  for  (8.6b) 
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Figure  ?b.    LW  sche're  for  (8.6c) 
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Fi-gur^  3c.     ULTl     scheme  for  (8.6c) 
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Figure  7a.  Hybrid  scheme  [(8.1)  +  (8.9)]  -  10  x  30  grid. 
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rig.  7b.  ULT2  scheme  -  10  x  30  grid. 
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